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We survey odd-even nuclear binding energy staggering using density functional theory with several 
treatments of the pairing interaction including the BCS, Hartree-Fock-Bogoliubov, and the Hartree- 
Fock-Bogoliubov with the Lipkin-Nogami approximation. We calculate the second difference of 
binding energies and compare with 443 measured neutron energy differences in isotope chains and 
418 measured proton energy differences in isotone chains. The particle-hole part of the energy 
functional is taken as the SLy4 Skyrme parametrization and the pairing part of the functional is 
based on a contact interaction with possible density dependence. An important feature of the data, 
reproduced by the theory, is the sharp gap quenching at magic numbers. With the strength of the 
interaction as a free parameter, the theory can reproduce the data to an rms accuracy of about 
0.25 MeV. This is slightly better than a single-parameter phenomenological description but slightly 
poorer than the usual two-parameter phenomenological form C/A". The following conclusions can 
be made about the performance of common parametrization of the pairing interaction: (i) there 
is a weak preference for a surface-peaked neutron-neutron pairing, which might be attributable to 
many-body effects; (ii) a larger strength is required in the proton pairing channel than in the neutron 
pairing channel; (iii) pairing strengths adjusted to the well-known spherical isotope chains are too 
weak to give a good overall fit to the mass differences. 



PACS numbers: 21.60.Jz, 21.30.Fe, 21.10.Dr 

I. INTRODUCTION 

The theory of nuclear masses or binding energies has 
attracted renewed interest with the advent of compu- 
tational resources sufficient to performed global calcu- 
lations based on self-consistent mean-field theory, also 
called density functional theory (DFT) 0,0,11. A long- 
term goal is an improved reliability for a theory that 
avoids ad hoc phenomenological parametrizations. One 
particular aspect of the nuclear binding problem is the 
ubiquitous phenomenon of odd-even staggering (OES) of 
binding energy. Since the early days of BCS theory [4] 
it has been largely attributed to BCS pairing, but there 
are in fact a number of mechanisms that can contribute 
d, i, 0, i, i, [13, EH, El- In this work we want to study 
the performance of BCS and its Hartree-Fock-Bogoliubov 
(HFB) extension to the global body of data, taking an 
energy density functional and pairing functionals that are 
in common use. In that way, we hope to provide a bench- 
mark to assess future improvements in the theory. Since 
we do not consider all mechanisms to generate the OES, 
our conclusions must be tentative. 

There are many DFT surveys that treat individual iso- 
tope chains, e.g., [l^, with the Z—50 isotope chain 
a favorite for calculation of pairing properties. We shall 
see, however, that it can be quite misleading to draw 
general conclusions without examining the whole body 
of OES data. Also, in much of the literature the OES 
was not obtained from differences of calculated binding 
energy but rather inferred from the average HFB gap 



parameters, as, e.g., in Ref. fl^. We also mention the 
global mass tables by the Brussels-Montreal collabora- 
tion [l^, • While this work achieves a good perfor- 
mance on binding energies, it deviates from the frame- 
work of DFT by adding phenomenological modifications 
to the theory. In particular, the pairing strength may 
depend on local densities but it is hard to justify an ex- 
plicit dependence on the number parity as is assumed in 
ref. dl. 

There are numerous measures of the OES in the liter- 
ature, inclu ding 3-point, 4-point, and 5-point difference 
formulas [3, [T9r[20l. [2ll [2^ . In this work we will use the 

3-point formula Ao as advocated in Ref. 6] and also 
used in Rcfs. [lO, EBl- For odd neutron number N, it is 
defined by the binding energy difference 

A(,3) (N) = i [B{N + 1) + B{N - 1) - 2B{N)] . (1) 

In the following, we shall call this quantity the neutron 
OES. Our survey will cover the proton OES as well. One 

(3) 

advantage of the Ao statistic is that it can be applied 
to more experimental data than the higher-order ones. 
Another advantage is that it suppresses the smooth con- 
tributions from the mean field to the gap. The other 3- 
point indicator, Ai'^'' (N) with A^-even, is less interesting 
for our purposes because it is more sensitive to single- 
particle energies. 

This paper is organized as follows. Section [Til outlines 
the theoretical DFT framework employed in this work. 
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In Sec. mil the selection of experimental data used in the 
survey is discussed. The results for selected spherical 
and deformed isotopic/isotonic chains are presented in 
Sec. lIVI while the global performance of our pairing mod- 
els is analyzed in Sec. El Finally, Sec. |Vl] contains the 
main conclusions and perspectives. 

II. METHODOLOGIES 

We carry out two independent surveys with the same 
Skyrme functional SLy4 [2^ in the particle-hole chan- 
nel. The pairing functional uses the zero-range density- 
dependent S interaction: 

V{ry)^Vo(^l-V^yir-v'). (2) 

Here Vq < is the pairing strength, p(r) is the isoscalar 
nucleonic density, and po^O.lGfm"'^. We have performed 
global calculations for ?7=0, 0.5, and 1, called volume, 
mixed, and surface, pairing, respectively. The volume 
pairing interaction acts primarily inside the nuclear vol- 
ume while the surface pairing generates pairing fields 
peaked around or outside the nuclear surface. As dis- 
cussed in Ref. [23|, different assumptions about the den- 
sity dependence can result in notable differences of pair- 
ing fields in neutron rich nuclei. 

The two surveys were carried out assuming two differ- 
ent theoretical frameworks for the pairing channel, the 
BCS and the Hartree-Fock-Bogohubov (HFB). The de- 
tails are described in the two subsections below. 



A. HF-BCS with ev8 

The HF-BCS extension of the nuclear DFT can be de- 
fined very concisely. The ordinary variables in the the- 
ory, namely the orbital wave functions expressed in 
some basis, are augmented by the BCS amplitudes Vi. 
Specifically, one defines the BCS Vi and Ui amplitudes 
for each orbital and calculates the ordinary DFT energy 
from its functional using the density matrix p(r,r') — 
t;?(/>*(r)(/)i(r'). To this is added the pairing energy 
functional, given by 

Epair = ^ VijUiViUjVj + ^ Vuvf (3) 
i^j i 

where Vij are the matrix elements of the pairing interac- 
tion. 

We use the code ev8 [Hi to carry out the HF+BCS 
computations. Ev8 solves the HF-t-BCS equations for 
Skyrme-type functionals via a discretization of the in- 
dividual wave- functions on a 3D Cartesian mesh and 
the imaginary time method. (In Ev8, the pairing func- 
tional ([3]) is approximated as jVijUiViUjVj.) The 
pairing interaction matrix elements are those of a density- 
dependent contact interaction Contact interactions 



can only be used in truncated orbitals spaces; the calcu- 
lations use the same truncation as is Ref. [1^ , namely an 
energy window of 10 MeV around the Fermi level. 

As we will see later, the OES only fluctuates about 
an average trend by ~0.3 MeV, putting a high demand 
on the accuracy and the nucleus-to-nucleus consistency 
of the self-consistent mean field calculations. The usual 
iteration procedure in EvS appears to be adequate to 
achieve accuracy at the 100 keV level in several hundred 
iterations at a fixed deformation. (Here accuracy means 
with respect to the fully converged minimum of the nu- 
merically implemented energy functional. This numer- 
ical functional may contain approximations that give a 
larger error with respect to the mathematically defined 
fimctional. In the case of Ev8, the lattice representation 
of the kinetic operator results in an error of the order of 
one MeV in heavy nuclei that varies very smoothly with 
A. Thus it largely cancels in the calculation of AJ, .) 
Finding the minima irrespective of deformation is less 
straightforward. We adopted the following protocol to 
determine them. We first build a table of DFT energies 
and orbital wave functions of the relevant even-even nu- 
clei, using the minimum energy deformations from the 
table calculated in Ref. [2^. The relative energies of 
spherical and deformed configurations are quite sensitive 
to the pairing interaction, so in the cases where the spher- 
ical configuration in that table has an energy close to the 
deformed minimum, the spherical was also tested. When 
it came out lower with the new pairing interaction, it re- 
placed the old entry in the new table. Next, the table 
was refined in an iterative way using only the deforma- 
tion information about neighboring even-even nuclei. If 
two neighboring nuclei have substantially different defor- 
mations in the table, each configuration must be tested 
in both nuclei. If taking the lower energy configuration 
results in a change, the process is repeated on the neigh- 
bors of the replaced nucleus. This is continued until there 
are no further changes in the even-even DFT solutions. 

Once the DFT table of even-even nucleus is finalized, 
the odd-A nuclei are calculated starting from the DFT 
solutions for the neighboring even-even nuclei. We per- 
formed the calculations using the so-called filling approx- 
imation for the odd particle [13, H^. The odd particle 
is assumed to occupy an orbital defined by its position 
in the list of orbitals ordered by single-particle energy. 
That orbital is blocked by setting v'^ = u'^ = 0.5 in the 
calculation of all ordinary densities, and omitting the 
orbital in the summation in Eq. ([3]). During the self- 
consistency iterations, the blocked orbital evolves along 
with the others, and thus may change character if the rel- 
ative ordering of the levels changes. Note that the filling 
approximation gives equal occupation numbers to both 
time-reversed partners, and therefore misses the effects 
of time-odd fields on the OES. 

Our protocol to find the most favorable orbital to block 
was to examine the five orbitals around the Fermi level 
of the neighboring even-even system. We also tested 
configurations generated from the DFT solution for the 
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even-even nucleus with one more nucleon than the target 
odd-A nucleus. Thus the total number of odd-A con- 
figurations tested is ten: five starting from the lighter 
even-even core and five starting from the heavier one. 

Since the objective is to determine the level of accu- 
racy that can be achieved, the calculations were carried 
for the nuclei in the data set for a number of values of 
the pairing strength Vq. The results below are reported 
for a value Vq close to that which minimizes the average 

(3) 

residual in the Ao data sets, taking neutron and protons 
independently. 



TABLE I: Effective strengths for the pairing interactions de- 
rived in the present study by means of Eq. (|4]).. Units of Vb 
are MeV-fm"^ . Also the Values of the fit parameter x are given 
in paratheses. 

theory density dependence V(^"{x) Vq''{x) 
'BCS volume 4651] 490.0 

mixed 700.0 755.0 

surface 1300.0 1462.0 

HFB mixed 318.1(0.41) 352.0(-0.18) 

HFB-f-LN mixed 300.5(0.18) 332.6(-0.44) 



B. HFB with hfbtho 

The HFB calculations were carried out with the ax- 
ial 2D HFB solver hfbtho [29,] that has recently been 
improved by implementing the modified Broyden mixing 
[301 to accelerate the convergence rate. 

The even-even nuclei are calculated first. An initial set 
of configurations is generated by performing constrained 
minimizations on a quadrupole deformation mesh. Typi- 
cally, there are 20 calculations having deformations in the 
range —0.5 < (3 < 0.5 with a mesh spacing A/3 = 0.05. 
Next, we turn off the constraint to find the local min- 
ima of the energy as a function of /3. When there are 
multiple minima, we select up to three for further pro- 
cessing, taking no more than one of oblate and prolate 
deformation, and also the spherical solution if it is a lo- 
cal minimum. The final step is to perform unconstrained 
minimizations on the selected configurations. The it- 
erative minimization is carried out until the maximum 
change of the matrix elements of the HFB matrix ele- 
ments falls below 0.0001 MeV. However, in one case the 
iteration converges to a limit cycle with energies oscillat- 
ing by 0.004 MeV. Since this is well below the accuracy 
needed here, we accepted the (lowest) calculated energy. 

The minimization for odd A?^ or Z is started from the 
candidate configurations produced at the second stage of 
the even-even calculations. As in the BCS, the odd nu- 
cleus is treated in the filling approximation, by blocking 
one of the orbitals. Here one has to specify which orbitals 
to block to generate the odd-nucleon configurations. The 
blocking candidates are determined by examining the 
HFB quasiparticle spectrum of the neighboring even-even 
nucleus with smaller number of nucleus. Tested are all 
one-quasiparticle configurations with quasi-particle ener- 
gies below the energy cutoff Eiqp_cut which is not smaller 
than 2 MeV for heavy nuclei and not bigger than 8 
MeV for very light systems. For most nuclei, we take 
-Eigp, cut =25/ vG4 MeV. As in the last step for even-even 
nuclei, unconstrained calculations are performed for all 
candidate configurations to find the absolute minimum 
energy. 

In the second variant of HFB calculations (HFB+LN) 
we performed an approximate particle number projection 
(before variation) using the Lipkin-Nogami (LN) method 
[sil . H^l ■ The practical implementation of the LN treat- 



ment follows Ref. [33| where the method was compared 
to the full particle number projection. 

In HFB and HFB-I-LN calculations we employed the 
orbital space extending to 20 major harmonic oscillator 
shells. For the pairing interaction, there is no lower en- 
ergy orbital cutoff and for the upper equivalent energy 
cutoff we adopted the commonly used value of 60 MeV 
[3^ . The calculations were first performed with a stan- 
dard pairing strength Vq*''^ adjusted to the average pair- 
ing gap in ^^°Sn according to the procedure of [33] ■ The 
values obtained are V^*"^ = -258.2 and -284.57 for the 
HFB and the HFB-f LN calculations. However, following 
our first survey, we found these strengths to be too small 
to make a good global fit to OES. We then increased the 
pairing strength by a factor of 1.2 and recalculated the 
mass table from scratch. Both sets of tables are available 
through the UNEDF SciDAC collaboration Our fit 
to the global systematics is then made using a linear fit 
of the data sets of the two mass tables: 

M{x) = xM{Vo'"^) + (1 - x)M{1.2 ■ Vq'"^). (4) 

The effective pairing strengths obtained in this way is 
given by 

Vq^^ {1.2-0.2x)V^"^. (5) 

The values of x and the derived pairing strengths are 
reported in Table [H 

C. Other methodological aspects 

In setting up the calculational protocols for the sur- 
vey, we of course scrutinized cases where the residuals 
between theory and experiment were large. The residuals 
are very sensitive to numerical inaccuracies, and their de- 
tailed analysis often brought to attention problems with 
calculated numbers due to, e.g., incorrectly assigned con- 
figurations or the lack of convergence in self-consistent 
iterations. We could then refine the protocols by mak- 
ing a broader screen or demanding higher precision to 
produce more accurate tables. From this standpoint we 
found Ao a very useful indicator. Also, the fact that the 
theory is variational is a tremendous help: any change in 
protocol that gives lower energies is necessarily an im- 
provement. 
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Apart from the numerics, it should be noted that 
the pairing gap is a very strong function of the pair- 
ing strength. For example, for the HFB calculations, we 
found that increasing the pairing strength by 20% from 
the value fitted to ^^"^Sn, the average neutron pairing gap 
increases by a factor of 2.3. This sensitivity is not surpris- 
ing. A typical nucleus in our set has some deformation 
and does not have a large single-particle degeneracy at 
the Fermi level. This implies that its BCS condensate 
is weak. Under these conditions, the BCS gap increases 
very quickly from a zero value at some finite strength of 
the interaction. In fact, we find that roughly 20-30% of 
the calculated OES contain a nucleus whose BCS or HFB 
condensate has collapsed. Even in infinite systems having 
no single-particle gap at the Fermi energy, the conden- 
sate A grows exponentially with interaction strength in 
the well-known BCS formula. 

An exact implementation of HFB requires the breaking 
of the time-reversal symmetry in the intrinsic frame of an 
odd-A nucleus. The resulting time-odd mean fields con- 
tribute to the binding energy of the odd-A system; hence, 

the OES aI^-* depends on whether these time-odd terms 
are included or not 0, H, [HI, [H, • In general, time- 
odd mean fields are known poorly and their effects differ 
from model to model. For instance, in the Skyrme DFT 
calculations of Ref. [sj , the time-odd polarizations sys- 
tematically shift the hole states down and particle states 
up in energy, while a different result has been obtained 
in the relativistic mean field approach 0, [1] . For the pur- 
pose of illustration. Fig. [1] shows the OES calculated in 
the HFB for isotonic chains including the effects of time- 
odd fields. A detailed analysis of resulting polarizations 
will be published elsewhere [38i] . It appears that the con- 
tribution is less that 100 keV in heavy elements such as 
the rare earths and actinides, but larger values have been 
found with other models, eg. [1, [l^l and in light nuclei. 
In any case, the density functional has not been fitted 
to time-odd properties, so no quantitative evaluation of 
their effect is possible. In this work, we employ the filling 
approximation, neglecting all time-odd interactions. 



III. EXPERIMENTAL DATA BASE 

The data base for our survey derives from the 2003 
atomic mass evaluation [40]. The accuracy requirements 
our purposes is of the order of 100 keV, so we only use 
the masses whose evaluated experimental errors were less 
than 200 keV. This gives 908 mass triplets for the neutron 
OES and 864 for protons. However, we made additional 
cuts to remove nuclei for which some physics is obviously 
missing from our BCS or HFB theory. First of all, in odd- 
odd nuclei there is an additional neutron-proton pairing 
effect. Its origin is not completely clear ^l|, but it is 
obviously beyond the scope of the pairing theory we use 
here. We therefore do not include binding energy triplets 
containing odd-odd nuclei. Another phenomenon in nu- 
clear binding that affects the OES is the so-called Wigner 
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FIG. 1: (Color online) Impact of time-odd terms on the pro- 
ton OES Ao^^ in the A''=90 and 100 isotonic chains calculated 
in the HFB theory with SLy4 functional and mixed pairing. 
Solid lines connect the values of proton OES including time- 
odd fields; dashed lines show the values of OES omitting the 
time-odd fields. Calculations were performed with the code 
HFODD of Ref. in a deformed HO basis containing 680 
deformed harmonic oscillator states. 



energy, an increased binding at N=Z. This might also be 
a neutron-proton pairing effect (see, e.g., Refs. [13, E^l)- 
We therefore eliminate triplets that contain N=Z nu- 
clei. More generally, mean-field approximations becomes 
doubtful when the number of particles is small. Some 
restriction of light nuclei is imposed by a cut on particle 
numbers, requiring that the neutron and proton numbers 
be greater than 8, corresponding to nuclei heavier than 
^^O. Finally, we only include nuclei with N>Z. There 
are only a few OES on the proton-rich side satisfying 
our other criteria, and keeping a fixed sign of the isospin 
will permit us to make some qualitative statement about 
the isospin dependence of the interaction. With all these 
cuts, there are left 443 triplets for the neutron OES and 
418 for the proton in our final data set. 

The sets of neutron and proton AJ, are plotted as a 
function of neutron and proton number in the top panels 
of Fig. [2l Lines connect the values of OES for the same 
number of nucleons of the opposite kind. It is common to 
plot OES as a function of A, but plotting it with respect 
to nucleons of the same kind makes clearer the origin of 
fluctuations. Probably the largest cause of fluctuation 
is the variation in single-particle level densities and the 
character of the level at the Fermi energy. This obviously 
depends strongly on the number of nucleons of the same 
kind, and this motivates the choice of abscissa variable. 
The single-particle level densities may also change with 
the different numbers of nucleons of the other kind, par- 
ticular if the additional nucleons causes a large change in 
deformation. Such effects should be visible in the varia- 
tion of OES at fixed value of the abscissa. To emphasize 
the like-nucleon fiuctuations, we also show as solid circles 
the average with respect to nucleons of the opposite kind. 



5 



One can see that the shell effects are large and there 
are large fluctuations on the scale of major shell spacings. 
For the neutron values of Ao , there are strong dips in 
the Z-averaged values at A'^=15, 29, 51, 83, and 125, i.e., 
in the vicinity of shell closures. We shall call this phe- 
nomenon gap quenching. Obviously, the OES is reduced 
when one of the three nuclei is at a magic number where 
the gap in single-particle energies is large. We will ex- 
amine the effect in more detail below, in presenting the 
theoretical OES. In Table [H we show the extreme OES 
cases- either the largest or the smallest in our experi- 
mental data sets. 



TABLE II: Nuclei from the data base adopted for our survey 
with the largest and smallest experimental values of A^'^' (in 
MeV). 





largest 




smallest 






{N,Z) 




{N,Z) 




neutrons 


(21,16) 


1.87 


(125,82) 


0.32 


protons 


(12,9) 


2.07 


(126,81) 


0.31 



Another observation that can be made about the neu- 
tron OES is that the variations with respect to Z are par- 
ticularly large in the regions iV=50-60 and 95-110. For 
the protons, the regions of strongest iV-dependence are 
A^=35-40 and 60-70. We will see that this is associated, 
at least in part, with changing deformations. 

The averages and variances for neutron and pro- 
ton OES aI^^ in the data set are 1.04±0.31 and 
0.96±0.27MeV, respectively. The lower average for the 
protons can be largely attributed to the Coulomb inter- 
action: in the liquid drop formula, the term acZ'^ / A^^^ 

f 31 

gives an average value of 0.11 McV for the proton Ao 
in the data set. There is some overall dependence of the 
OES on mass number A which may be seen by visual 
inspection of Fig. O For more discussion of the global 
mass dependence of the OES, we refer the reader to 
Refs. |6j, tlSi, l22] . The smooth yl-dependence seems rather 
weak compared to the local nucleus-to-nucleus fluctu- 
ations caused by shell effects, but parameterizing it in 
some way can give improved fits. For example, the phe- 
nomenological parametrization [l9l . [i^ l 



with c=4.66MeV (4.31 MeV) and a=0.31 gives an rms 
residual of 0.25 MeV on the neutron (proton) data set. 
The global trends given by Eq. ([6]) are shown as the 
dashed lines in Fig. O For the sake of the plot, we av- 
eraged over nucleon number of the opposite kind just as 
was done to produce experimental averages. 

An important question about the global systematics 
is whether there is an isospin dependence of the pair- 
ing interaction. It is clear from Fig. [5] that there can be 
strong interaction between the pairing of one kind on the 



numbers of the other kind, but as mentioned it could be 
due to other effects such as shape changes. The isospin 
dependence has been examined in Ref. (2^ but no sig- 
nificant effect has been found (see also Refs. [1^, [4^ H^l 
for more discussion of isovector trends). On the other 
hand, a recent study [1^ of the OES of nuclear masses 
for isotopic chains between the proton shell closures at 
Z—50 and Z—82, including nuclei with extreme isospins, 
has claimed a significant isospin dependence of pairing. 
We will discuss in more detail the possible evidence for 
an isovector dependence of the interaction in Sec. IV CI 



IV. RESULTS: LOCAL COMPARISONS 

We begin our comparison between theory and exper- 
iment with two spherical semi-magic isotope chains, Sn 
and Pb. The results of the calculations are shown in 
Fig. [31 For all three treatments of pairing, the trends 
of predicted OES are consistent with the data, concern- 
ing both global and local variations. In particular, the- 
ory reproduces the flatness in the Sn isotopes up to the 
quenched gap at A^ = 83, and the downsloping trend 
from the light Pb isotopes up to the quenched gap at 
^"^Pb. In the Sn isotope chain, there is a small dip at 
A^ = 65 which might be attributed to a neutron subshell 
closure at A^ = 64. In any case, the theories all predict 
a shallow local minimum. When confronted with experi- 
ment, HFB appears to do better slightly than HF-I-BCS. 
As mentioned earlier, the strength of the pairing interac- 
tion was fit to the overall systematics, giving a somewhat 
too high an average OES in both spherical chains. The 
performance of the theories for long isotonic and isotopic 
chains that include deformed nuclei is illustrated in Fig. 3] 
and [5l Fig. 2] shows the neutron OES in the Dy isotopes 
(Z=66) and proton OES in the A^=98 isotones. Like in 
the case of semi-magic nuclei, the agreement with experi- 
ment is good, in particular for well deformed nuclei where 
the mean field changes smoothly with particle number. 
The effect of changing deformation is illustrated by the 
region from A ~ 160 to ~ 190 which starts deformed and 
becomes spherical as Z is increased from 66 to 82. The 
neutron values of OES for A^ = 102 covering this transi- 
tion region are shown in Fig. [5] as a function of Z. The 
squares show the experimental OES, which increase from 
about 0.6-0.8 MeV for the lower Z nuclei and goes up to 
'--^ 1.3 MeV for the singly magic Z = 82 isotope. The 
circles show the corresponding calculated HFB values of 
OES. The trend is very similar, but the theoretical rise to 
Z=82 is sharper and higher than is seen experimentally. 
Very likely, the increase in single-particle level density 
going from deformed to spherical nuclei is responsible for 
increasing trend in the OES. 

The fact that our calculations overestimate OES in 
spherical nuclei may be partly attributed to the parti- 
cle number fluctuations. The pairing gap exponentially 
depends on the inverse of the single-particle level density 
at the Fermi level, which is large in spherical open-shell 
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FIG. 2: OES as a function of nucleon number. From top to bottom, panels show Aj, for experiment, HFB, HFB+LN, and BCS 
treatments, respectively. Circles show values obtained by averaging over nucleon number of the opposite isospin to that of the 
OES. The calculation used the SLy4 Skyrme energy functional and a pairing interaction with the mixed density dependence. 
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FIG. 3: (Color online) Calculated (HF+BCS, HFB, and 
HFB+LN) and experimental values of aI^^ for neutrons in 
the Sn (top) and Pb (bottom) semi-magic isotopic chains. 
In all calculations, the pairing interaction was taken in the 
mixed pairing form (?7=0.5) with strength Vo (or x in Eq. Q) 
adjusted to global systematics. 



FIG. 4; (Color online) Calculated and experimental values of 
Ao"^^ for neutrons (top, Z—66) and protons (bottom, A^=98) 
in rare earth nuclei, including strongly deformed systems. 
Filled circles: experiment; squares: BCS; open circles HFB; 
triangles: HFB+LN. 



nuclei due to the 2j+l degeneracy (the limit of pair- 
ing rotation). In deformed systems, the level density is 
reduced due to the Jahn- Teller effect [H, [s^l and this 
gives rise to the overall pairing reduction; hence, a transi- 
tion towards the transitional pairing regime, in which the 
particle number fluctuations are more important. Since 
the original pairing strength Vq^'^ was adjusted to the 
global data set containing far more non-spherical nuclei 
than semi-magic systems, the resulting deformation bias 
results in too strong pairing correlations predicted for 
spherical nuclei, as seen in Fig. [31 This can be partly 
cured by considering particle- number fluctuations. In- 
deed, as seen in Figs. [3] and S] the LN procedure shghtly 
improves agreement with experiment for spherical nuclei 
while still reproducing data for deformed chains. 
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FIG. 5: Experimental (filled circles) and HFB (open circles) 
neutron OES for A'' — 102 as a function of Z. The chain starts 
in the well-deformed lanthanides and ends next to the singly 
magic ^**''Pb. 



RESULTS: GLOBAL 



Overview and shell effects 



The lower three sets of panels in Fig. [2] show distri- 
butions of neutron and proton OES for our three theo- 
retical treatments: HFB, HFB-I-LN and BCS. To make 
the overall trends with respect to shell filling more ap- 



parent, we also show the values of the OES obtained by 
averaging over nucleons of the opposite isospin. Qualita- 
tively, the three methodologies give rather similar results. 
In all cases, the trends of the predicted gaps are consis- 
tent with the data. The strong neutron gap quenching 
seen in the experimental OES at the numbers TV = 83 
and N = 125 is reproduced in all three theoretical treat- 
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TABLE III: Characteristics of nuclei with quenched gaps in 
the HFB calculations. The listed nucleus is the one with 
smallest OES at the given gap. The calculated deformation (3 
is given in the third column. The last columns give the quasi- 
particle quantum numbers angular momentum and parity 
for spherical nuclei and azimuthal angular momentum K and 
parity K-n for deformed nuclei. 



N gap 


Nucleus 


13 


r or K'' 


15 


^^Ne 


0.00 


1/2+ 


29 


52rpj 


0.08 


1/2- 


47-51 


»^Kr 


-0.09 


5/2+ 


83 


""Gd 


-0.03 


7/2" 


125 




0.00 


1/2- 


Z gap 


15 


39p 


0.21 


1/2+ 


29 


61 Cu 


0.10 


3/2- 


47-49 


"^Ag 


-0.22 


1/2+ 


69 


"^Tm 


0.33 


7/2- 


81 


203 rpj 


0.01 


1/2+ 



ments. The gap quenching is obviously dependent on the 
presence of shell closures, but the fact that it does not 
invariably occur on both sides of the magic numbers indi- 
cates that particular orbital properties must play a role. 
Since the HFB calculations were performed in an axially 
symmetric basis, we can examine the quantum numbers 
of the blocked orbital. In Table III we show the quasi- 
particle orbital characteristics for odd nuclei exhibiting 
quenched gaps. The large quenching at N — 125 can 
be understood as a spherical shell effect associated with 
the pi/2 shell at the Fermi level. A j = 1/2 shell will 
have reduced pairing for two reasons. First, there is no 
degeneracy within the shell to correlate pairs, so all of 
the pairing has to come from off-diagonal interactions to 
other shell. Second, these couplings are reduced because 
the spatial overlap of high-/ and low-Z orbitals is poor. 
Similar considerations apply to A^=15 and Z = 81, where 
the relevant spherical shell is 

While the location of the gap quenching is well re- 
produced by theory, the magnitude of the effect is often 
exaggerated. Most notably, all the experimental values 
of OES are positive, the theoretical ones at {N^ Z) = 
(20 — 24, 15) even have a negative sign. 

Comparing the different treatments of pairing, we see 
that the fluctuations seen are at the same positions in 
BCS and HFB-fLN as in the HFB. However, the am- 
plitudes of the fluctuations seem somewhat larger in the 
BCS treatment but somewhat smaller in HFB+LN. It 
is not clear why the BCS should emphasize the fluctu- 
ations, but the fact that they are damped in HFB+LN 
is not surprising. In both BCS and HFB the static pair- 
ing sometimes collapses in a significant fraction of nuclei, 
as mentioned earlier. The pairing never collapses in the 
HFB-I-LN treatment, so the OES should be smoother as 
one passes into a region of weak pairing. 

Large variations with proton number are found around 
N^7Q and in the region iV=100-120. As mentioned ear- 
lier, the latter region is a transition region between spher- 



ical and deformed ground states, and that will affect the 
OES. The HFB theory has a spike as well as quenched- 
gap behavior at Z = 81 in the proton gap. The origin of 
the spike appears to be the coexistence of spherical and 
deformed configurations in light isotopes of Z—81 and 82 
nuclides. At the phase transition point, there can be a 
large static polarization contribution to the OES. 

To see where theory performs best, and where impor- 
tant physics is missing, in Fig. [6] we show the isospin- 
averaged residuals for our HFB model. As expected, the 
best agreement is obtained for well deformed rare earths 
and actinides whose properties vary smoothly with parti- 
cle number. The largest deviations are seen around shell 
closures and in the regions of shape coexistence (^^--^90 
for neutrons and ^=110 and 190 for protons) where dy- 
namic shape fluctuations are known to strongly impact 
masses (26l. IsTjl. 



0.6 - 



0.4 



0.2 



> 

u 



Neutron variances (exp - HFB) 



20 40 60 80 100 120 140 160 180 200 220 240 

Mass number A 



0.8 - 




80 100 120 140 160 180 200 220 240 

Mass number A 

FIG. 6: (Color online) Neutron (top) and proton (bottom) 
average variances a between HFB and experimental Ai^' as 
a function of ^. At each A, the variance a was obtained by 
averaging the residuals over all isobars available. 

To analyze more quantitatively the improvement in ac- 
curacy for deformed nuclei, we have examined the rms 
residuals for the neutron OES separating the data into 
bins by the calculated values of the deformation (3. Fig- 
ure [7] shows the averaged HFB residuals. As expected, 
the residuals gradually decrease with deformation. The 
transitional/coexisting nuclei with weakly-oblate shapes 
show the largest deviations from the data. 
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Neutron variances (HFB-exp) 



-0.1 0.1 0.2 0.3 

quadrupole deformation 

FIG. 7: Neutron average variances a between HFB and ex- 
perimental Ao'^' as a function of the calculated quadrupole 
deformation (52- The deformation range -0.2</32<0.4 was di- 
vided into bins with A/32=0.1, and each bin the variance a 
was obtained by averaging the residuals over all nuclei avail- 
able. 



The fluctuations in OES are obviously suppressed 
when one averages over nucleons of the opposite isospin. 
With that averaging, the comparison between theory and 
experiment looks much better. Fig. [8] shows the theory- 
experiment comparison for the HFB methodology. Be- 
sides seeing the shell effects discussed above, one also sees 
better how the theory performs with respect to the A- 
dependence of the pairing. The theoretical proton OES 
has an overall A-dependence that seems to accord well 
with the experimental trend. For the neutron OES, how- 
ever, the theory is flatter than the experimental trend. 

We have carried out the survey with different assump- 
tions about the density dependence of the pairing inter- 
action to see sensitivity of the A-dependence to that char- 
acteristic. The results of the averaged neutron OES for 
volume and surface pairing are shown in Fig. [5] The 
effect is very small, except for the light nuclei. In view 
of the other possible contributions to the staggering, we 
do not believe that one can reliably extract the density 
dependence of the effective pairing interaction strength 
from the observed ^-dependence. We discuss below an- 
other mechanism that could simulate the observed trend 
in A, namely an isospin dependence of the effective pair- 
ing interaction. 



B. Performance statistics 



In Table IrVl we report the rms residuals for the OES in 
the various treatments, fitting the strength of the pairing 
interaction separately for neutrons and protons. In the 
case of our HFB and HFB-I-LN models, we carried out 
separate optimizations for neutrons and protons with re- 
spect to the X parameter in Eq. ^ . The effective pairing 
strength obtained in this way is given by Eq. ([5]). A more 
proper procedure would be to make a two-dimensional 
optimization based on recalculated HFB mass tables as- 
suming different strengths of proton and neutron pairing. 
Our experience with HF-I-BCS model, however, is that 
the neutron pairing strength does not significantly affect 

(3) 

the proton A), and vice versa, at the level of changes 
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FIG. 8: (Color online) Comparison between calculated (HFB 
with mixed pairing) and experimental aI^' values. Top: Z- 
averaged values for neutrons. Bottom: TV-averaged values for 
protons. 



considered in this work. Therefore, for the purpose of a 
global survey, a simplified treatment has been adopted. 

From Table [TVl we see that all three treatments of the 
pairing can achieve the performance of zeroth order de- 
scription as a constant (or phenomenological) gap, but 
only the Lipkin-Nogami, shown on the next-to-last line, 
does significantly better. In the HFB-I-LN, pairing corre- 
lations are always present. This is particularly important 
for odd-v4 nuclei, where the standard blocking procedure 
often gives rise to the unphysical pairing collapse, arti- 
ficially affecting the OES and producing an exaggerated 
fluctuation. However, one should be cautious in using the 
HFB+LN. While for open-shell systems using it gives a 
good agreement with those of the HFB with the full par- 
ticle number projection before variation, the method is 
inaccurate for closed-shell systems [s^. Consequently, 
it is safest to use HFB+LN only away from the magic 
numbers. In addition, the numerical procedure to find 
the solution lacks stability when there is a large gap at 
the Fermi level. Nevertheless, we obtained converged so- 
lutions for 440 out the 443 neutron triplets and 411 out 
the 418 proton triplets with our HFB+LN implementa- 
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N 

FIG. 9: Averaged BCS neutron OES with the volume (squares 
connected by dotted Une) and surface (circles connected by 
solid line) pairing interactions of Eq. ([5}. The dot-dashed 
curve is the phenomenological fit to the data using the func- 
tional form ((6)1. 



tion. The numbers in Table HVl are for those data sets. If 
we restrict the data set further to omit the magic num- 
bers 28,50,82, and 126, the rms residual of the neutron 
OES is hardly affected, changing from 0.23 to 0.22. 



TABLE IV: RMS residuals of A^''' obtained in various mod- 
els. All energies are in MeV. The last column shows the ra- 
tio of proton and neutron effective pairing strengths obtained 
through the optimization procedure. The mass predictions of 
the IIFB-14 model 17] were taken from 52]. 



Theory 


pairing 


residual 


residual 








neutrons 


protons 




Constant 




0.31 


0.27 




c/A" 




0.24 


0.22 




HF+BCS 


volume 


0.31 


0.38 


1.05 


HF-HBCS 


mixed 


0.30 


0.36 


1.08 


HF+BCS 


surface 


0.27 


0.35 


1.12 


HFB 


mixed 


0.27 


0.33 


1.11 


HFB+LN 


mixed 


0.23 


0.28 


1.11 


HFB-14 




0.46 


0.44 


1.10 



One of the basic questions about nuclear pairing is the 
role of induced interactions in the effective pairing in- 
teraction [H, [13, m, [5^ . Indirect information about 
this can in principle be obtained by exhibiting the den- 
sity dependence and the isospin dependence of the ef- 
fective interaction. It is therefore of interest to examine 
interactions including a density dependence to see the 
sensitivity. The rms residual for the neutron OES with 
volume, mixed, and surface pairing in HF+BCS theory 
are shown in Table IIVI There is a slight favoring of the 
surface interaction, but we deem that the difference in 
the residuals ( 10 %) is too slight to be significant. The 
weak sensitivity to the density dependence confirms the 
results of other studies [l^, [53| ■ 



C. Isospin dependence 

A possible isospin dependence of the effective pair- 
ing strength has been much discussed in the literature 
[2a . m, 1^, 113, m, HI- The nuclear interaction may 
be assumed to conserve isospin at a fundamental level 
but the coupling to core excitations can be different for 
neutron and protons when the core has a neutron ex- 
cess. Another isospin-dependent contribution to pairing 
comes from the Coulomb interaction. Indeed, inclusion 
of the Coulomb has been found to substantially suppress 
the pairing interaction energy [6^ and the pairing gaps 
[sst . In the last column of the Table IIVI we report the 
ratio of neutron and proton interaction strengths we ex- 
tract from our fits to Ao^-* . The effective proton strength, 

(3) 

needed to reproduce experimental Ao , is larger than the 
neutron strength. If the Coulomb were included explic- 
itly, we would expect that the needed nuclear interaction 
would be even larger for the protons. Since the under- 
lying strong interaction is isoscalar to a good approxi- 
mation, we believe that our inferred isospin effect must 
arise from induced three-body interactions involving the 
neutron excess. We note in passing that a number of 
mass table fits by the Goriely group arrive at pairing 
strengths larger for protons than neutrons. An exam- 
ple is the HFB-14 model [l^, shown in the last line of 
Table ITVl However, since different interactions are used 
for even and odd nuclei in HFB-14, the results are not 
directly comparable. 

As another way to test the data for an isospin- 
dependent pairing interaction, we separate the nuclei into 
two subsets according to neutron excess, and compare the 
average residual OES. To define the subsets, we first di- 
vide the nuclei into five A-bins. For each bin we make 
a cut at some value of J = (TV — Z) /A to have roughly 
equal numbers for the two sets, which we designate "low 
isospin" and "high isospin" . In that way the effect of any 

A-dependence in the OES will be reduced. The binning 

(3) 

for proton and neutron values of Ao has to be done sep- 
arate to get balanced sets. The average values of OES 
for the two sets are reported in Table |Vl The empir- 



TABLE V: Average A)f' (in MeV) calculated in HFB+LN 
sorted by neutron excess. See text for details. 



Data set 


low isospin 


high isospin 


difference 


neutrons exp 


1.13 


0.94 


-0.19 


HFB+LN 


1.05 


1.02 


-0.03 


protons exp 


1.05 


0.88 


-0.17 


HFB+LN 


0.99 


0.93 


-0.06 



ical OES is lower for higher neutron excesses for both 
protons and neutrons. The calculated Ao for neutrons 
are nearly equal, while the calculated Ao for protons 
do show a difference but much smaller than observed. 
For both cases, the differences would require weakening 
the pairing interaction for nuclei with larger neutron ex- 
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cesses. An isospin dependence would have the opposite 
sign for proton and neutron gaps. 

As a final plausibility check on whether the different 
strengths could be attributable to some isoscalar three- 
body interaction, we computed the OES for nuclei on 
the opposite side of the N = Z line (recall that our fitted 
data set is restricted to > Z). This comprises 5 neu- 
tron OES and 6 proton OES, excluding as before cases 
involving N — Z nuclei. Taking the pairing strengths 
from our global fit, we find that the calculated average 
neutron OES is low (by 0.7 MeV) while the calculated 
average proton OES is too high (by 0.24 MeV). This is 
precisely the expected direction of the error if the dif- 
ference in the effective strengths depends on the sign of 
N~Z, as would be required by an overall isoscalar energy 
functional 

Thus, the main evidence for an isospin dependence in 
the present theory is the need for different strengths for 
the overall fits to neutron and proton data sets. This 
results supports the recent attempts [H, to directly 
parametrize the pairing functional in terms of isovector 
densities. 



VI. PERSPECTIVE 

The present study demonstrates that the current state 
of the art in the nuclearDFT permits calculation of OES 
to accuracy of the order of 0.25 MeV rms. This is not a 
trivial outcome in two respective. First, the binding ener- 
gies involved range up to nearly four orders of magnitude 
larger, so there is a high demand on computational pre- 
cision in carry out the DFT. Second, the pairing gap is a 
highly sensitive function of the mean field properties such 
as level density, and so the theory needs have an accurate 
treatment of the single-particle properties. In this work 
we have only considered the SLy4 functional which has 
known deficiencies. Clearly further exploration of DFT 
functionals is warranted, perhaps including ones having 
different (isoscalar and isovector) effective masses 

We found no large differences between the BCS and 
HFB treatments. This is not unexpected; it is only near 
the drip lines that the HFB with its better treatment 
of spacial variations of the anomalous density is needed. 
On the other hand, we believe that the improvement we 
found for the HFB-I-LN treatment is significant, showing 
that a number-conserving treatment of the pairing cor- 
relations is needed. Many of the nuclei, particularly the 
odd ones, are on the edge of collapse of BCS pairing, and 
for these a number-conserving treatment is essential to 
calculate the pairing correlation energy. Unfortunately, 
the LN treatment of number violation is not reliable near 
closed shells. We would therefore advocate in future us- 
ing other treatments of number violation, perhaps HFB 
with variation-after-projection ^33.] , or mapping onto an 
exactly soluble pairing Hamiltonian [6lj . 

The global data on OES shows a weak A-dependence 



that is certainly not reproduced with a pairing interac- 
tion that is density-independent or has only a mild depen- 
dence on density. In the calculations with the BCS the- 
ory, the sensitivity to density dependence was explored, 
and it was found that there were only small changes in 
the overall performance. It should also be noted that the 
mean-field contributions to the OES are highly depen- 
dent on A. Thus, firm conclusions about the origin of 
the A dependence must await surveys based on theory 
that avoids the filling approximation. 

A very interesting question, related to the density de- 
pendence, is where there is an effective isospin depen- 
dence of the strength of the pairing interaction. The 
question cannot be addressed with confidence by exam- 
ining individual isotope or isotone chains, because the 
other species of nucleon can affect the effective pairing 
Hamiltonian, particularly the single-particle spectrum. 
However, from the global survey we find what seems to 
be a robust result, that the effective pairing strength for 
protons OES is about 10% larger than for neutron OES. 
The calculation does not take into account the Coulomb 
interaction in the pairing channel, but naively that would 
be expected to decrease the effect strength, not increase 
it. The other possibility to explain the difference is an 
induced isospin dependence. 

It is perhaps disappointing that the overall perfor- 
mance of the theory is only slightly better than the naive 
one-parameter phenomenology attributing the staggering 
to a constant BCS gap. The two-parameter phenomeno- 
logical function c/A" does slightly better than the theory 
overall. But that form has no justification and the shell 
effects that are faithfully reproduced by the theory are 
missing. So we conclude that the rms residuals between 
theory and experiment do not tell the whole story. In 
any case, the promising possibility to surpass the per- 
formance of the present phenomenology is to continue 
the DFT studies with better functionals, including mean 
field contributions and number-conserving treatments of 
pairing. 
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